function impedance_from_N(CONSTS, plot_data)
        
    do_interpolation = false;

    phi = 0;
    m_max = 20; % 50
    N_low_order = 11;
    N_high_order = 13;
    num_points = 5;
    refine_multiplier = 10;
    
    Z_color = 'b';
    Z_odn_color = 'r';
    
    N_vec = logspace(N_low_order, N_high_order, num_points)';
    N_refined_vec = logspace(N_low_order, N_high_order, refine_multiplier*num_points)';
    Z_vec = zeros(size(N_vec));
    Z_odn_pl_vec = zeros(size(N_vec));
    
    for i = 1:size(N_vec,1)
        N = N_vec(i);

        LOCAL_CONSTS = CONSTS;
        LOCAL_CONSTS.N = N;
        LOCAL_CONSTS = update_CONSTS(LOCAL_CONSTS);

        I_full = func_I_full(phi, m_max, LOCAL_CONSTS);
        Z_vec(i) = 1 ./ I_full;

        I_odn_pl = sum_I_odn_pl(phi, LOCAL_CONSTS);
        Z_odn_pl_vec(i) = 1 ./ I_odn_pl;
    end
    
    coef_sgs_to_si = 9e+11;
    Z_sgs_vec = abs(coef_sgs_to_si*Z_vec);
    Z_odn_sgs_vec = abs(coef_sgs_to_si*Z_odn_pl_vec);
    
    Z_sgs_refined_vec = spline(N_vec, Z_sgs_vec, N_refined_vec);
    Z_odn_sgs_refined_vec = spline(N_vec, Z_odn_sgs_vec, N_refined_vec);

    if (plot_data)
        figure;
        
        if(do_interpolation)
            semilogx( ...
                N_vec, Z_sgs_vec, [Z_color, '.'], ...
                N_refined_vec, Z_sgs_refined_vec, [Z_color, '-'], ...
                N_vec, Z_odn_sgs_vec, [Z_odn_color, 'x'], ...
                N_refined_vec, Z_odn_sgs_refined_vec, [Z_odn_color, '-'] ...
                );
        else
            semilogx( ...
                N_vec, Z_sgs_vec, [Z_color, '.-'], ...
                N_vec, Z_odn_sgs_vec, [Z_odn_color, 'x-'] ...
                );
        end  
        
        grid; 
        title('Abs of Z(N)'); 
        xlabel('N, cm^{-3}'); ylabel('|Z(N)|, Om');
        %legend('|Z(N)|','|Z(N)| for homogen. plasma', ... % legend is broken
        %    'location', 'NorthEast');
    end

end